Análise de Variância Multivariada (MANOVA)

ANOVA \(\times\) MANOVA


A análise multivariada de variância (MANOVA) é uma generalização da análise de variância univariada (ANOVA)

Enquanto a ANOVA se concentra no estudo das diferenças entre as médias populacionais dos grupos em uma única variável dependente, a MANOVA examina essas diferenças em duas ou mais variáveis dependentes simultaneamente.

ANOVA \(\times\) MANOVA

Tanto a ANOVA quanto a MANOVA só fornecem informações sobre a existência de diferenças estatisticamente significativas entre as médias dos grupos em um conjunto de variáveis dependentes, necessitando de outros testes para obtermos resultados mais confiáveis e precisos.

A maioria das pesquisas não estão interessadas em avaliar as diferenças entre médias para uma única variável dependente, mas sim para um conjunto de variáveis dependentes.

Se fizermos várias ANOVAs, presumimos que não existe uma estrutura de correlação entre as variáveis dependentes, ou que uma tal estrutura não é de interesse.

ANOVA \(\times\) MANOVA

A MANOVA presume existência de correlação significativa entre as variáveis dependentes.

Neste contexto, se realizarmos várias ANOVAs

ERRO TIPO I: Rejeitar \(H_0\) quando esta é verdadeira!

ANOVA: Pequena revisão…

Modelo ANOVA a um fator univariado


Considere g amostras aleatórias independentes \(X_{j1}, X_{j2}, \cdots, X_{jn_j}\), \(j = 1, \cdots, g\) de tamanhos \(n_j\) de distribuições \(N(\mu_j,\sigma^2)\).

As hipóteses de interesse são:

\[H_0: \mu_1 = \mu_2 = \cdots = \mu_g\] \[H_a: \text{pelo menos uma média é diferente das demais}\]

Modelo ANOVA a um fator univariado

O modelo ANOVA é definido por

\[X_{ji} = \mu_j + \epsilon_{ji}\]

\[X_{ji} = \underbrace{\mu_j}_{\text{média da j-ésima amostra}} + \overbrace{\epsilon_{ji}}^{\text{erro aleatório i-ésima observação na j-ésima amostra}}\]

\(i = 1, \cdots, n_j\) e \(j = 1, \cdots, g\)

Modelo ANOVA a um fator univariado

Em geral, adota-se a seguinte reparametrização:

\[\mu_j = \mu + \tau_j\]

\[\mu_j = \underbrace{\mu}_{\text{média geral}} + \overbrace{\tau_{j}}^{\text{efeito do j-ésimo tratamento}}\]

e, assim, as hipóteses equivalentes são:

\[H_0: \tau_1 = \tau_2 = \cdots = \tau_g = 0\] \[H_a: \text{pelo menos um é não-nulo}\]

Modelo ANOVA a um fator univariado

O modelo ANOVA é redefinido por

\[X_{ji} = \mu + \tau_j + \epsilon_{ji}\]

\(i = 1, \cdots, n_j\) e \(j = 1, \cdots, g\)

Suposição: \(\epsilon_{ji} \stackrel{\text{ind}}{\sim} N(0,\sigma^2)\)

Modelo ANOVA a um fator univariado

No modelo de efeito fixo, temos:

\[\mu = \displaystyle{\frac{\displaystyle{\sum_{j=i}^g n_j \mu_j}}{\displaystyle{\sum_{j=1}^g n_j}}} = \displaystyle{\frac{\displaystyle{\sum_{j=i}^g n_j (\mu + \tau_j)}}{\displaystyle{\sum_{j=1}^g n_j}}} = \displaystyle{\mu + \frac{\displaystyle{\sum_{j=i}^g n_j \tau_j}}{\displaystyle{\sum_{j=1}^g n_j}}} \Rightarrow \displaystyle{\sum_{j=1}^g n_j\tau_j = 0}\]

Modelo ANOVA a um fator univariado

O procedimento básico adotado é a decomposição da soma de quadrados totais:

\[\begin{eqnarray} \underbrace{\displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(x_{ji} - \bar{x}_{\ldotp \ldotp})^2}_{\text{SQ Total}} &=& \displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(x_{ji} - \bar{x}_{j \ldotp} + \bar{x}_{j \ldotp} - \bar{x}_{\ldotp \ldotp})^2 = \nonumber \\ &=& \underbrace{\displaystyle{\sum_{j=1}^g} n_j (\bar{x}_{j \ldotp} - \bar{x}_{\ldotp \ldotp})^2}_{\text{SQ Tratamentos}} + \underbrace{\displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(x_{ji} - \bar{x}_{j \ldotp})^2}_{\text{SQ Resíduos}} \nonumber \end{eqnarray}\]

Modelo ANOVA a um fator univariado

As estimativas de mínimos quadrados de \(\mu\) e \(\tau_j\) são dadas por

\[\widehat{\mu} = \bar{x}_{\ldotp \ldotp} = \displaystyle{\frac{1}{n}}\displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}} x_{ji},\text{ com }n = \displaystyle{\sum_{j=1}^g n_j}\]

\[\widehat{\tau}_j = \bar{x}_{j \ldotp} - \bar{x}_{\ldotp \ldotp}, \text{ com } \bar{x}_{j \ldotp} = \displaystyle{\frac{1}{n_j}} \displaystyle{ \sum_{i=1}^{n_j} x_{ji}}\]

Um estimador não viciado para \(\sigma^2\) é dado por \(\widehat{\sigma}^2 = \displaystyle{\frac{\text{SQ Resíduos}}{n - g}}\), com \(n = \displaystyle{\sum_{j=1}^g n_j}\).

Modelo ANOVA a um fator univariado

Tabela ANOVA a um fator em g níveis


Fonte de Variação G.L. SQ QM F
Tratamento \(g - 1\) \(\displaystyle{\sum_{j=1}^g} n_j (\bar{x}_{j \ldotp} - \bar{x}_{\ldotp \ldotp})^2\) \(\displaystyle{\frac{SQ Trat}{g - 1}}\) \(\displaystyle{\frac{QM Trat}{QM Res}}\)
Resíduos \(n - g\) \(\displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(x_{ji} - \bar{x}_{j \ldotp})^2\) \(\displaystyle{\frac{SQ Res}{n - g}}\)
Total \(n - 1\) \(\displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(x_{ji} - \bar{x}_{\ldotp \ldotp})^2\)

Modelo ANOVA a um fator univariado


Se \(H_0\) é verdadeira, \(F \sim F_{g - 1, n-g}\). Logo, rejeitamos \(H_0\) ao nível de significância \(\alpha\) se \(F > F_{g - 1, n-g}(\alpha)\).

Observe que isto é equivalente a dizer que o teste rejeita \(H_0\) para valores grandes da razão \(\displaystyle{\frac{SQ Trat}{SQ Res}}\).

Modelo ANOVA a um fator univariado


Equivalentemente, valores grandes de \(1 + \displaystyle{\frac{SQ Trat}{SQ Res}}\), ou de forma análoga, valores pequenos da recíproca, dada por

\[\displaystyle{\frac{SQ Res}{SQ Trat + SQ Res}}\]

Adiante, veremos que uma estatística de teste multivariada, com suas devidas adaptações, tem forma similar à esta acima.

MANOVA: Uma generalização da ANOVA

Modelo MANOVA a um fator


Considere as seguintes amostras aleatórias coletadas de \(g \geqslant 2\) populações:

Amostra 1: \(\mathbf{x}_{11}, \mathbf{x}_{12}, \cdots, \mathbf{x}_{1n_1}\) da população \(\pi_1\)

Amostra 2: \(\mathbf{x}_{21}, \mathbf{x}_{22}, \cdots, \mathbf{x}_{2n_2}\) da população \(\pi_2\)

\[\vdots\]

Amostra g: \(\mathbf{x}_{g1}, \mathbf{x}_{g2}, \cdots, \mathbf{x}_{gn_g}\) da população \(\pi_g\)

Modelo MANOVA a um fator


  • Suposições
  1. \(\mathbf{x}_{j1}, \mathbf{x}_{j2}, \cdots, \mathbf{x}_{jn_j}\) é uma amostra de tamanho \(n_j\) de uma população com média \(\mathbf{\mu}_j\), \(j = 1, 2, \cdots, g\);
  1. Todas as populações têm matriz de covariâncias comum \(\mathbf{\Sigma}\);
  1. Cada população é normal p-variada;
  1. As amostras são independentes.

Modelo MANOVA a um fator

O modelo MANOVA é definido por

\[\mathbf{x}_{ji} = \mathbf{\mu} + \mathbf{\tau}_j + \mathbf{\epsilon}_{ji}\]

\(i = 1, \cdots, n_j\) e \(j = 1, \cdots, g\)

  • Suposição: \(\mathbf{\epsilon}_{ji} \stackrel{\text{ind}}{\sim} N_p(\mathbf{0},\mathbf{\Sigma})\)

  • Restrição: \(\displaystyle{\sum_{j=1}^g n_j \mathbf{\tau}_j = \mathbf{0}}\)

Modelo MANOVA a um fator


Soma de Quadrados Total:

\[SQ Tot = \displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(\mathbf{x}_{ji} - \bar{\mathbf{x}}_{\ldotp \ldotp})(\mathbf{x}_{ji} - \bar{\mathbf{X}}_{\ldotp \ldotp})^t\]

Modelo MANOVA a um fator

  • Observações:
  1. \(\bar{\mathbf{x}}_{\ldotp \ldotp}\) é o vetor de médias amostral geral;
  1. Agora, SQ Tot é uma matriz \(p \times p\);
  1. Os elementos de sua diagonal principal correspondem às somas de quadrados totais para cada variável separada;
  1. Fora da diagonal, temos as somas dos produtos cruzados dos desvios, caracterizando a estrutura de dependência das p variáveis estudadas.

Modelo MANOVA a um fator

Decomposição da Soma de Quadrados Total:

\[\underbrace{\displaystyle{\sum_{j=1}^g} n_j (\bar{\mathbf{x}}_{j \ldotp} - \bar{\mathbf{x}}_{\ldotp \ldotp}) (\bar{\mathbf{x}}_{j \ldotp} - \bar{\mathbf{x}}_{\ldotp \ldotp})^t}_{\text{SQ Tratamentos}} + \underbrace{\displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(\mathbf{x}_{ji} - \bar{\mathbf{x}}_{j \ldotp})(\mathbf{x}_{ji} - \bar{\mathbf{x}}_{j \ldotp})^t}_{\text{SQ Resíduos}}\]

Observação: Os graus de liberdade de todos os casos são os mesmos do caso univariado.

  • Hipóteses:

\[H_0: \mathbf{\tau}_1 = \mathbf{\tau}_2 = \cdots = \mathbf{\tau}_g = \mathbf{0}\] \[H_a: \text{pelo menos um é não-nulo}\]

Modelo MANOVA a um fator

Tabela MANOVA a um fator em g níveis


Fonte de Variação G.L. Matriz de SQ
Tratamento \(g - 1\) \(H = \displaystyle{\sum_{j=1}^g} n_j (\bar{\mathbf{x}}_{j \ldotp} - \bar{\mathbf{x}}_{\ldotp \ldotp}) (\bar{\mathbf{x}}_{j \ldotp} - \bar{\mathbf{x}}_{\ldotp \ldotp})^t\)
Resíduos \(n - g\) \(E = \displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(\mathbf{x}_{ji} - \bar{\mathbf{x}}_{j \ldotp}) (\mathbf{x}_{ji} - \bar{\mathbf{x}}_{j \ldotp})^t\)
Total \(n - 1\) \(T = \displaystyle{\sum_{j=1}^g} \displaystyle{\sum_{i=1}^{n_j}}(\mathbf{x}_{ji} - \bar{\mathbf{x}}_{\ldotp \ldotp}) (\mathbf{x}_{ji} - \bar{\mathbf{x}}_{\ldotp \ldotp})^t\)

Modelo MANOVA a um fator

  • Observações:
  1. \(T = H + E\)
  1. As diagonais principais de cada uma dessas matrizes de somas de quadrados e produtos cruzados fornecerão as somas de quadrados da ANOVA para cada variável separadamente.

Uma possível estatística de teste envolve variâncias generalizadas. Seja,

\[\Lambda = \displaystyle{\frac{|E|}{|H + E|}}\]

Rejeitamos \(H_0\) se \(\Lambda\) for um valor pequeno.

Modelo MANOVA a um fator


  • A estatística \(\Lambda\) foi originalmente proposta por Wilks e corresponde a uma forma equivalente ao teste F da hipótese de ausência de efeito de tratamento do caso univariado.
  • A estatística lambda de Wilks tem a vantagem de ser conveniente e estar correlacionada ao teste da razão de verossimilhança.
  • distribuição exata de \(\Lambda\) pode ser derivada para os casos especiais apresentados a seguir:

Modelo MANOVA a um fator

Casos especiais

\(p = 1\) \(g \geqslant 2\) \(\displaystyle{\left(\frac{n - g}{g - 1}\right)} \displaystyle{\left(\frac{1 - \Lambda}{\Lambda}\right)} \sim F_{g - 1,n - g}\)
\(p = 2\) \(g \geqslant 2\) \(\displaystyle{\left(\frac{n - g - 1}{g - 1}\right)} \displaystyle{\left(\frac{1 - \sqrt{\Lambda}}{\sqrt{\Lambda}}\right)} \sim F_{2(g - 1),2(n - g - 1)}\)
\(p \geqslant 1\) \(g = 2\) \(\displaystyle{\left(\frac{n - p - 1}{p}\right)} \displaystyle{\left(\frac{1 - \Lambda}{\Lambda}\right)} \sim F_{p, n - p- 1}\)
\(p \geqslant 1\) \(g = 3\) \(\displaystyle{\left(\frac{n - p - 2}{p}\right)} \displaystyle{\left(\frac{1 - \sqrt{\Lambda}}{\sqrt{\Lambda}}\right)} \sim F_{2p,2(n - p - 2)}\)
  • Observação: Para outros casos e tamanhos amostrais grandes, uma modificação de \(\Lambda\) devido da Bartlett pode ser usada para testar \(H_0\).

Modelo MANOVA a um fator


Bartlett mostrou que se \(H_0\) é verdadeira e \(n\) é grande, então,

\[-\left(n - 1 - \displaystyle{\frac{p + g}{2}}\right) \ln(\Lambda) \dot{\sim} \chi_{p(g - 1)}^2\]

Consequentemente, se \(n\) é grande, rejeitamos \(H_0\) ao nível de significância \(\alpha\), se

\[-\left(n - 1 - \displaystyle{\frac{p + g}{2}}\right) \ln(\Lambda) > \chi_{p(g - 1)}^2(\alpha)\]

Modelo MANOVA a um fator

Outros testes multivariados

  • Traço de Hotelling-Lawley:

\[U = tr(HE^{-1})\]

Sob \(H_0\), temos que

\[\displaystyle{\frac{2(sn + 1)}{s^2(2m + s + 1)}U \hspace{0.2cm} \dot{\sim} \hspace{0.2cm} F_{s(2m + s + 1),2(sn + 1)}}\]

Modelo MANOVA a um fator

Outros testes multivariados

  • Traço de Pillai:

\[V = tr[H(H + E)^{-1}]\]

Sob \(H_0\), temos que

\[\left(\displaystyle{\frac{V}{s - V}} \right) \left(\displaystyle{\frac{2n + s + 1}{2m + s + 1}}\right) \hspace{0.2cm} \dot{\sim} \hspace{0.2cm} F_{s(2m + s + 1),s(2n + s + 1)}\]

Modelo MANOVA a um fator

Outros testes multivariados

  • Raiz máxima de Roy:

\[\Theta = \lambda_1\]

sendo \(\lambda_1\) o maior autovalor da matriz \(HE^{-1}\).

Sob \(H_0\), temos que

\[\displaystyle{\frac{\Theta(\nu - d + q)}{d}} \hspace{0.2cm} \dot{\sim} \hspace{0.2cm} F_{d,\nu - d + q}\]

Modelo MANOVA a um fator

Outros testes multivariados


Observações

\(\nu\) = G.L. Resíduo \(d = \max(p,q)\) \(s = \min(p,q)\)
\(n = \displaystyle{\frac{\nu - p - 1}{2}}\) \(m = \displaystyle{\frac{|p - q| - 1}{2}}\) \(q\) = G.L. Tratamento
\(p\) = número de variáveis

Modelo MANOVA a um fator

Exemplo - Análise dos dados de uma casa de repouso para idosos:

O departamento de saúde e serviços sociais de determinada cidade reembolsa casas de repouso para idosos no estado por serviços oferecidos.

O departamento desenvolve um conjunto de fórmulas para as taxas com facilidade, baseado em fatores tais como níveis de cuidados, taxa média de salários e taxa média de salários no estado.

As casas de repouso podem ser classificadas com relação à propriedade: privadas, sem fins lucrativos e públicas e também com relação à certificação: especializadas em enfermagem, unidades de cuidados intermediários ou uma combinação dos dois.

Modelo MANOVA a um fator

O objetivo foi verificar os efeitos de propriedade e certificação (ou ambos) sobre os custos. Quatro custos, calculados por paciente/dia, medidos em horas/paciente, foram usados:

  • \(Y_1\): Custo de mão de obra de enfermagem;

  • \(Y_2\): Custo de nutricionista;

  • \(Y_3\): Custo de trabalho de manutenção e funcionamento;

  • \(Y_4\): Custo de limpeza e lavanderia.

Um total de \(N = 516\) observações sobre os \(p = 4\) custos foram separadas pelo tipo de propriedade e estão disponíveis no arquivo Exemplo Manova.dat.

Modelo MANOVA a um fator

pacman::p_load(
  tidyverse,
  rstatix,
  ggpubr
)
dados = read.table("https://raw.githubusercontent.com/tiagomartin/est022/refs/heads/main/dados/Exemplo_Manova.dat", sep = ' ',header = TRUE) %>% mutate(tipo = as.factor(tipo))
dados %>% str()
'data.frame':   516 obs. of  5 variables:
 $ y1  : num  2.14 2.19 2.12 2.75 2.59 ...
 $ y2  : num  0.5 0.543 0.58 0.367 0.492 0.362 0.494 0.379 0.621 0.389 ...
 $ y3  : num  0.122 0.07 0.135 0.075 0.014 0.113 0.122 0.117 0.05 0.06 ...
 $ y4  : num  0.246 0.459 0.501 0.466 0.258 0.443 0.302 0.282 0.38 0.411 ...
 $ tipo: Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...

Modelo MANOVA a um fator

Com fins ilustrativos assumiremos normalidade e igualdade de variâncias.

# H0: mu1 = mu2 = mu3

modelo = manova(cbind(y1, y2, y3, y4) ~ tipo, data = dados)

## Teste lambda de Wilks
summary(modelo, test = "Wilks")
           Df   Wilks approx F num Df den Df    Pr(>F)    
tipo        2 0.79914   15.126      8   1020 < 2.2e-16 ***
Residuals 513                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
H = summary(modelo)$SS$tipo

E = summary(modelo)$SS$Residuals

Lam = det(E)/det(H+E)
Lam
[1] 0.7991428

Modelo MANOVA a um fator

## Teste traço de Pillai 
summary(modelo, test = "Pillai")
           Df  Pillai approx F num Df den Df    Pr(>F)    
tipo        2 0.20726   14.769      8   1022 < 2.2e-16 ***
Residuals 513                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Teste raiz máxima de Roy
summary(modelo,test="Roy")
           Df     Roy approx F num Df den Df    Pr(>F)    
tipo        2 0.20406   26.069      4    511 < 2.2e-16 ***
Residuals 513                                             
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Teste traço de Hotelling-Lawley
summary(modelo, test = "Hotelling-Lawley")
           Df Hotelling-Lawley approx F num Df den Df    Pr(>F)    
tipo        2          0.24333   15.482      8   1018 < 2.2e-16 ***
Residuals 513                                                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Intervalos simultâneos de confiança para efeitos de tratamento

Quando \(H_0\) é rejeitada, os efeitos que levaram à rejeição são de interesse.

Para efeito de comparações pareadas, a abordagem de Bonferroni pode ser usada para construir intervalos simultâneos de confiança para os componentes dos vetores diferença \(\mathbf{\tau}_j - \mathbf{\tau}_k\).

Seja \(\tau_{ji}\) o i-ésimo componente de \(\mathbf{\tau}_j\). Como \(\mathbf{\tau}_j\) é estimado por \(\widehat{\mathbf{\tau}}_j = \bar{\mathbf{x}}_{j \ldotp} - \bar{\mathbf{x}}_{\ldotp \ldotp}\), temos que o i-ésimo componente de \(\mathbf{\tau}_j\) é estimado por \(\bar{{x}}_{ji} - \bar{{x}}_{\ldotp i}\).

Intervalos simultâneos de confiança para efeitos de tratamento

Assim, \(\tau_{ji} - \tau_{ki}\) é estimado por

\[\widehat{\tau}_{ji} - \widehat{\tau}_{ki} = \bar{x}_{ji} - \bar{x}_{ki}\]

que é uma diferença entre médias de duas amostras independentes.

O teste t para duas amostras independentes é válido para um nível de significância apropriado.

Intervalos simultâneos de confiança para efeitos de tratamento

Temos ainda que,

\[\text{Var}(\widehat{\tau}_{ji} - \widehat{\tau}_{ki}) = \left(\displaystyle{\frac{1}{n_j}} + \displaystyle{\frac{1}{n_k}}\right) \sigma_{ii}\]

com \(\sigma_{ii}\) o i-ésimo elemento da diagonal de \(\Sigma\). Como sugerido, \(\sigma_{ii}\) é estimado pelo i-ésimo elemento da diagonal da matriz \(\mathbf{E}\) dividido pelos respectivos graus de liberdade, a saber, \(\widehat{\sigma}_{ii} = \displaystyle{\frac{1}{n - g}}e_{ii}\).

Intervalos simultâneos de confiança para efeitos de tratamento

Resta identificar o número de intervalos. Como são \(g\) tratamentos para \(p\) medidas, segue que ao todo teremos \(m = \displaystyle{\dfrac{pg(g-1)}{2}}\) intervalos.

Logo, com confiança de pelo menos \((1 - \alpha)\), os intervalos de confiança simultâneos para \(\tau_{ji} - \tau_{ki}\) são dados por:

\[(\bar{x}_{ji} - \bar{x}_{ki}) \pm t_{n-g} \left(\displaystyle{\frac{\alpha}{pg(g - 1)}}\right) \sqrt{\left(\displaystyle{\frac{1}{n_j}} + \displaystyle{\frac{1}{n_k}}\right)\displaystyle{\frac{e_{ii}}{n - g}} }\]

em que \(e_{ii}\) é o i-ésimo elemento da diagonal da matriz \(\mathbf{E}\).

Intervalos simultâneos de confiança para efeitos de tratamento

Exemplo - Análise dos dados de uma casa de repouso para idosos

summary.aov(modelo)
 Response y1 :
             Df  Sum Sq Mean Sq F value    Pr(>F)    
tipo          2   7.473  3.7363  9.9001 6.044e-05 ***
Residuals   513 193.608  0.3774                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response y2 :
             Df Sum Sq Mean Sq F value    Pr(>F)    
tipo          2 1.0430 0.52149  33.268 2.601e-14 ***
Residuals   513 8.0415 0.01568                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response y3 :
             Df  Sum Sq  Mean Sq F value    Pr(>F)    
tipo          2 0.13858 0.069290   25.48 2.816e-11 ***
Residuals   513 1.39502 0.002719                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 Response y4 :
             Df Sum Sq  Mean Sq F value    Pr(>F)    
tipo          2 0.4048 0.202389  15.299 3.517e-07 ***
Residuals   513 6.7863 0.013229                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Intervalos simultâneos de confiança para efeitos de tratamento

dados %>% emmeans_test(y1 ~ tipo, p.adjust.method = "bonferroni", detailed = TRUE)
# A tibble: 3 × 14
  term  .y.   group1 group2 null.value estimate     se    df conf.low conf.high
* <chr> <chr> <chr>  <chr>       <dbl>    <dbl>  <dbl> <dbl>    <dbl>     <dbl>
1 tipo  y1    1      2               0   -0.172 0.0642   513   -0.298   -0.0458
2 tipo  y1    1      3               0   -0.295 0.0701   513   -0.433   -0.157 
3 tipo  y1    2      3               0   -0.123 0.0791   513   -0.279    0.0324
# ℹ 4 more variables: statistic <dbl>, p <dbl>, p.adj <dbl>, p.adj.signif <chr>
dados %>% emmeans_test(y2 ~ tipo, p.adjust.method = "bonferroni", detailed = TRUE)
# A tibble: 3 × 14
  term  .y.   group1 group2 null.value estimate     se    df conf.low conf.high
* <chr> <chr> <chr>  <chr>       <dbl>    <dbl>  <dbl> <dbl>    <dbl>     <dbl>
1 tipo  y2    1      2               0  -0.106  0.0131   513  -0.132    -0.0807
2 tipo  y2    1      3               0  -0.0447 0.0143   513  -0.0728   -0.0167
3 tipo  y2    2      3               0   0.0617 0.0161   513   0.0300    0.0934
# ℹ 4 more variables: statistic <dbl>, p <dbl>, p.adj <dbl>, p.adj.signif <chr>

Intervalos simultâneos de confiança para efeitos de tratamento

dados %>% emmeans_test(y3 ~ tipo, p.adjust.method = "bonferroni", detailed = TRUE)
# A tibble: 3 × 14
  term  .y.   group1 group2 null.value estimate      se    df conf.low conf.high
* <chr> <chr> <chr>  <chr>       <dbl>    <dbl>   <dbl> <dbl>    <dbl>     <dbl>
1 tipo  y3    1      2               0 -0.0287  0.00545   513  -0.0394  -0.0180 
2 tipo  y3    1      3               0 -0.0370  0.00595   513  -0.0487  -0.0253 
3 tipo  y3    2      3               0 -0.00835 0.00672   513  -0.0215   0.00484
# ℹ 4 more variables: statistic <dbl>, p <dbl>, p.adj <dbl>, p.adj.signif <chr>
dados %>% emmeans_test(y4 ~ tipo, p.adjust.method = "bonferroni", detailed = TRUE)
# A tibble: 3 × 14
  term  .y.   group1 group2 null.value estimate     se    df conf.low conf.high
* <chr> <chr> <chr>  <chr>       <dbl>    <dbl>  <dbl> <dbl>    <dbl>     <dbl>
1 tipo  y4    1      2               0  -0.0633 0.0120   513 -0.0869    -0.0397
2 tipo  y4    1      3               0  -0.0426 0.0131   513 -0.0684    -0.0168
3 tipo  y4    2      3               0   0.0207 0.0148   513 -0.00840    0.0498
# ℹ 4 more variables: statistic <dbl>, p <dbl>, p.adj <dbl>, p.adj.signif <chr>

Intervalos simultâneos de confiança para efeitos de tratamento

Intervalos de confiança simultâneos com correção de bonferroni de 95% para os grupos 1 e 2

Grupo variáveis dif médias lim inf lim sup
1 -0,172 -0,298 -0,046
1,2 2 -0,106 -0,132 -0,081
3 -0,029 -0,049 -0,025
4 -0,063 -0,087 -0,039

Intervalos simultâneos de confiança para efeitos de tratamento

Intervalos de confiança simultâneos com correção de bonferroni de 95% para os grupos 1 e 3

Grupo variáveis dif médias lim inf lim sup
1 -0,295 -0,433 -0,157
1,3 2 -0,045 -0,073 -0,017
3 -0,037 -0,049 -0,025
4 -0,043 -0,068 -0,017

Intervalos simultâneos de confiança para efeitos de tratamento

Intervalos de confiança simultâneos com correção de bonferroni de 95% para os grupos 2 e 3

Grupo variáveis dif médias lim inf lim sup
1 -0,123 -0,279 0,032
2,3 2 0,062 0,030 0,093
3 -0,008 -0,021 0,005
4 0,021 -0,008 0,049

Intervalos simultâneos de confiança para efeitos de tratamento

Observações

  • Codificação: 1 - Estabelecimento privado; 2 - Estabelecimento sem fins lucrativos; 3 - Estabelecimento público
  • Foram destacados os intervalos que excluíram o zero.
  • Em relação aos estabelecimentos privados, podemos notar que estes apresentam custos totais mais baixos em relação aos outros tipos de estabelecimento.
  • Considerando os estabelecimentos sem fins lucrativos e os públicos, notamos que não existem diferenças entre os custos, exceto para custos com nutricionistas, que parece ser maior em entidades sem fins lucrativos.

Visualizando graficamente

pwc <- dados %>%
  gather(key = "variables", value = "value", y1, y2,y3,y4) %>%
  group_by(variables) %>%
  games_howell_test(value ~ tipo)
pwc
# A tibble: 12 × 9
   variables .y.   group1 group2 estimate conf.low conf.high    p.adj
 * <chr>     <chr> <chr>  <chr>     <dbl>    <dbl>     <dbl>    <dbl>
 1 y1        value 1      2       0.172   -0.00643    0.351  6.2 e- 2
 2 y1        value 1      3       0.295    0.161      0.429  1.39e- 6
 3 y1        value 2      3       0.123   -0.0696     0.316  2.9 e- 1
 4 y2        value 1      2       0.106    0.0711     0.142  6.30e-11
 5 y2        value 1      3       0.0447   0.0124     0.0771 4   e- 3
 6 y2        value 2      3      -0.0617  -0.105     -0.0187 2   e- 3
 7 y3        value 1      2       0.0287   0.0136     0.0437 3.74e- 5
 8 y3        value 1      3       0.0370   0.0218     0.0522 1.63e- 7
 9 y3        value 2      3       0.00835 -0.0119     0.0286 5.96e- 1
10 y4        value 1      2       0.0633   0.0325     0.0941 6.94e- 6
11 y4        value 1      3       0.0426   0.0128     0.0724 3   e- 3
12 y4        value 2      3      -0.0207  -0.0579     0.0165 3.89e- 1
# ℹ 1 more variable: p.adj.signif <chr>

Visualizando graficamente

pwc <- pwc %>% add_xy_position(x = "tipo")
test.label <- create_test_label(
  description = "MANOVA", statistic.text = quote(italic("F")),
  statistic = 15.13, p= "<0.0001", parameter = "8,1020",
  type = "expression", detailed = F
)
ggboxplot(
  dados, x = "tipo", y = c("y1", "y2", "y3", "y4"), 
  merge = TRUE, palette = "jco"
) + 
  stat_pvalue_manual(
    pwc, hide.ns = TRUE, tip.length = 0, 
    step.increase = 0.1, step.group.by = "variables",
    color = "variables"
  ) +
  labs(
    subtitle = test.label,
    caption = get_pwc_label(pwc, type = "expression")
  )

Visualizando graficamente